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ABSTRACT 

High-resolution observations of the inner regions of barred disk galaxies have revealed 
many asymmetrical, small-scale central features, some of which are best described as 
secondary bars. Because orbital timescales in the galaxy center are short, secondary 
bars are likely to be dynamically decoupled from the main kiloparsec-scale bars. Here, 
we show that regular orbits exist in such doubly-barred potentials and that they can 
support the bars in their motion. We find orbits in which particles remain on loops: 
closed curves which return to their original positions after two bars have come back 
to the same relative orientation. Stars trapped around stable loops could form the 
building blocks for a long-lived, doubly-barred galaxy. Using the loop representation, 
we can find which orbits support the bars in their motion, and what are the constraints 
on the sizes and shapes of self-consistent double bars. In particular, it appears that a 
long-lived secondary bar may exist only when an Inner Lindblad Resonance is present 
in the primary bar, and that it would not extend beyond this resonance. 

Key words: celestial mechanics, stellar dynamics — galaxies: kinematics and dy- 
namics — galaxies: nuclei — galaxies: spiral — galaxies: structure 



1 INTRODUCTION 



Recent high-resolution imaging, from space and from the 
ground, shows hitherto-unsuspected small-scale structures 
at the centers of many galaxies. In a number of barred galax- 
ies, isophotal twists have been seen within the central few 
hundred parsecs. Small central bars, and occasional triply 
barred systems, are found: see Friedli (1996) for a review. 
The inner bars appear to be oriented randomly with respect 
to the larger bars (Buta & Crocker 1993), as expected if 
they are dynamically distinct subsystems. High-resolution 
CO maps of some disk galaxies show inner bars of molecu- 
lar gas {e.g. Devereux et al. 1992, Benedict et al. 1996), but 
the presence of inner bars in infrared images (see e.g. Friedli 
et al. 1996; Jungwiert, Combes & Axon 1997; Mulchaey, 
Regan & Kundu 1997) suggests they contain old stars, and 
do not consist purely of young stars and gas. A high fre- 
quency of double bars is seen in a multicolor imaging survey 
of early-type barred galaxies recently completed by Erwin & 
Sparke (1999). The survey, performed with the WIYN tele- 
scope, and supplemented with archival HST images, covers 
a reasonably complete sample of 38 nearby {z <2000 km/s), 
bright, barred SO and Sa galaxies in the field. At least ~20% 
of these galaxies appear to possess secondary bars. In this 
paper, we will call the larger bar the big, main, primary or 
outer bar. The smaller bar will be referred to as the small, 
secondary or inner bar. Like Friedli (1996) , we prefer to avoid 
the term 'nuclear bar', since a small bar on scales of tens or 



hundreds of parsecs can exist even in the absence of a main 
kiloparsec-scale bar. 

Within about 100 pc of a galactic center, orbital times 
are at least an order of magnitude less than those at a few 
kiloparsecs; thus a dynamically decoupled inner bar is likely 
to rotate faster than the outer structure. The entire pat- 
tern is then not steady in any reference frame. Such de- 
coupled secondary bars have been seen to form in numeri- 
cal simulations involving gravitating particles together with 
dissipative 'gas clouds' (Friedli & Martinet 1993, Heller & 
Shlosman 1994, Combes 1994, Shaw et al. 1995). The only 
systematic study of the orbital dynamics in a double-bar is 
that of Pfenniger & Norman (1990); they explored a weakly 
dissipative form of the equations of motion for a particle in 
one double-bar potential, finding a limit cycle in the plane 
of the bars, and spheroidal strange attractors. Orbits in the 
doubly barred potential do not have a conserved integral of 
motion, and in principle they might all be chaotic, exploring 
large regions of phase space. If the orbits are mostly chaotic, 
it is unlikely that such a system could be self-consistent, so 
that the average density of all the stars on their orbits in the 
time-varying potential adds up exactly to what is needed to 
give rise to the potential in which they move. 

How can potentials including two independently ro- 
tating bars maintain themselves as gravitating systems? 
What are the conditions under which a gravitationally self- 
consistent double-bar structure could exist? We approach 
these questions by considering particle orbits in models that 
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include two rigid bars rotating at two constant, incommen- 
surable pattern speeds. Such models cannot be fully self- 
consistent; Louis & Gerhard (1988) have shown that two 
bars which make up a self-consistent system must distort 
slightly as they rotate through each other. Nevertheless, we 
are looking for density distributions that are most likely to 
be close to a self-consistent model, i.e. for those which sup- 
port orbital families capable of hosting sets of particles that 
together recreate the assumed time-dependent density dis- 
tribution. In our investigation we use the concept of the 
loop introduced by Maciejewski & Sparke (1997). We ex- 
plore constraints on double-bars that could maintain a likely 
self-consistent structure. 

The concept of the loop is presented in §2. In §3, we re- 
view the methods that we used to find and to follow loops nu- 
merically in a doubly barred potential. In §4, we describe our 
models for the density distribution in doubly barred galax- 
ies, constructed on the basis of parameters derived from ob- 
servations. Various examples of loop families are given in 
§5. In §6, we closely examine loops supporting a likely self- 
consistent model of a double bar. Limitations imposed by 
orbital structure on self-consistent doubly barred systems, 
and their implications for gas flows, are discussed in §7. We 
summarize our results in §8. Analytical techniques that al- 
low us to approximate the loops supporting a doubly barred 
density distribution, and to model the gas streamlines, are 
presented in Appendix A. In Appendix B, we derive the for- 
mulae for a potential of a prolate Ferrers bar. 



2 THE CONCEPT OF THE LOOP 

A gravitational potential consisting of two concentric bisym- 
metric bars that are independently rotating is not stationary 
in any reference frame. Nevertheless, if the bars rotate in the 
same plane around their centers with angular velocities 57 b , 
fls, the potential pulsates with a frequency tOp = 2{fls — ^B) 
in the frame rotating with one of the bars (note that the fac- 
tor 2 comes from bisymmetry). 

In the stationary potential of a single bar, particles on 
a closed periodic orbit move along it, always staying on the 
same curve. Stable periodic orbits form "backbones" of a 
steady potential: nearby orbits are trapped around them. 
Orbits in a double-barred galaxy will generally not be closed 
in any uniformly rotating frame of reference, since particles 
there undergo two forcing actions with non-commensurable 
frequencies. We want to extend the definition of an orbit, in 
order to find closed curves which can similarly serve as back- 
bones of a non-steady, doubly barred system. We postulate 
that these are curves which, when populated with particles 
moving in a doubly barred potential, will return to their 
original positions every time the bars come back to the same 
relative orientation. We call these curves loops — they are a 
generalization of the closed periodic orbit in a single bar (see 
Maciejewski & Sparke 1997). A particle that begins its orbit 
from a position along a given loop returns to another point 
on the same loop after the bars have realigned. Loops change 
their shape as the bars rotate through each other. Particles 
trapped around loops that stay aligned with the bars in their 
motion could build up a long-lived, self-consistent, doubly- 
barred galaxy, in the same way as particles trapped near 
closed periodic orbits in a singly barred potential. 



Although here applied to two rigid bars rotating with 
constant pattern speeds, this method works for any potential 
pulsating with a constant period. In particular, the two bars 
do not have to be rigid, or to rotate with constant pattern 
speeds. 



3 METHODS 

3.1 Formal statement of the loop problem 

We consider a phase space S, whose points are the coordi- 
nates and momenta of particles at a given time. Then we 
define the mapping ViS), which assigns to each point in S 
the coordinates of the corresponding particle after integrat- 
ing its motion for one relative period Tp — 27r/ajp of the 
bars. This mapping P(5) can be applied repeatedly, leading 
to consecutive iterated maps. For each point A of the phase 
space S, the V"{A) is called the n^^ iterate of this point. 
The set of all iterates of a given point is called the orbit of 
this point under the mapping V{S) (see e.g. Lichtenberg & 
Lieberman 1992). Thus defined, the orbits divide the phase 
space into non-overlapping invariant sets: no point can be 
on two distinct orbits. 

A Hamiltonian system in A'' dimensions with a period- 
ically varying potential can be written as an autonomous 
system (one with a time-independent potential) in A'^ -I- 1 
dimensions (see e.g. Lichtenberg & Lieberman 1992, Louis 
& Gerhard 1988). If we look only at orbits in the symmetry 
plane in which the bars rotate, our 2-D potential changing 
periodically with a fixed frequency is equivalent to a 3-D 
autonomous Hamiltonian system. Regular orbits in this 3- 
D system are those that behave as if they have 3 integrals 
of motion: they are confined to a 3-D hypersurface of the 
6-D extended phase space of the autonomous system. Our 
iterative mapping suppresses the dimension corresponding 
to time, so the 3-D hypersurface of a regular orbit is repre- 
sented by a 2-D set of points. Thus if a point A lies on a 
regular orbit in the double-bar system, its orbit under ViS) 
is a 2-dimensional hypersurface in the 4-D space. Chaotic 
orbits will occupy 3 or 4 dimensions. If a point in the phase 
space S lies on an orbit which closes after A'^ times the rel- 
ative period, its orbit under the mapping V{S) will consist 
of a 0-dimensional set of N points. 

What about a point A £ S, whose iterates form a 1- 
D curve? This point is not on a periodic orbit (since these 
form a 0-dimensional set under 'P{S)), but it corresponds 
to a particle moving on a special kind of a regular orbit 
(those are only confined to 2-dimensional hypersurfaces). 
Every time the two bars return to the same relative position, 
the iterates 'P"(A) always lie on this 1-D curve: we have a 
loop as defined in §2. A particle that starts from any position 
on the loop will in general return to another place on the 
same loop every time the bars align. Since no point can be 
on two distinct orbits, and since in general the iterates of 
one point on the loop fill the entire loop, every point yl(x,v) 
which lies on a loop contains all the information about this 
loop. 

In Figure |l|a, we show how consecutive iterates oi A € S 
generate a loop in a doubly barred galaxy. Figure |i[d shows 
200, 400 and 4000 iterates 7^(5) for two points: the first hes 
near a stable loop, the second one does not. For the near-loop 
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Figure 1. a) A loop in position space, generated by iterates of 
one of its points: 5 first iterates (left), 200 iterates (middle) and 
4000 iterates (right), b) 200, 400 and 4000 iterates (left to right) 
in position space of a point that lies near a stable loop (top) and 
one that does not (bottom). Points are plotted when both bars 
are horizontal. The corresponding diagrams in the velocity space 
have the same appearance. 

solution, the points are confined to a ring around the stable 
loop. Iterates of a typical point far from the loop quickly dis- 
perse, covering much of the phase space. Thus a stable loop 
can be found by adjusting the coordinates of point A £ S, 
until its orbit under P(5) converges onto a 1-D closed curve. 
Unstable loops can only be found in this way if we have a 
very good guess at the starting conditions. Our method re- 
quires following a single particle for many bar rotations, and 
any particle starting close to an unstable loop will be lost 
by then. But since only stable loops can support the doubly 
barred systems, this is not a serious disadvantage. Another 
way to look at stability is to compare the ring widths for e.g. 
100 and 200 iterations. They should be similar for a stable 
loop. For an unstable loop, a nearby particle moves away 
exponentially with time, increasing the width. 

3.2 Numerical methods of finding the loop 

To find iterates of point A G S,we follow a particle in the xy 
plane in a given potential $(x, t), in which two bars rotate 
independently about the 2;-axis. We solve the equation of 
motion in a reference frame rotating with pattern speed fis, 
which is equal to the angular speed of one of the bars. In 
the xy plane of the barred disk, fts -L r, and one can write 
the equation of motion (see e.g. formula 3-82 in Binney & 
Tremaine 1987) as 

r = -V$ - 2(^3 X r) + I^bI^ r. (1) 

In Cartesian coordinates, it can be decomposed into its x 
and y parts: 

X = -— + 20.BV^^\x, (2) 
ox 

9$ 2 

y = IVlBX + VlBy. (3) 

oy 

We integrated the equations above using a variable-order, 
variable-step Adams integrating subroutine for first order 



differential equations. This predictor-corrector method is 
well suited to our problem, for which we expect multiply 
periodic solutions. Under the forcing action of two bars, a 
particle can be given a kinetic energy large enough to es- 
cape. Any particle departing from the potential center to 
distances a few times larger than the extent of the outer bar 
was considered to be lost, and we looked only for loops in 
and around the bars. 

We expect symmetric loops to be most important in 
symmetric bars. Therefore we restricted our search to the 
loops which are symmetric with respect to both x- and y- 
axes when the bars are aligned on the a;-axis: we reflected 
the successive iterates ^"(A) about the two axes so that 
they occupy the < <^ < 7r/2 region only. For symmetric 
loops the t/-velocity on the y-axis is zero, and one can ex- 
plore various initial i-velocities for a given y-axis cross- 
ing coordinate y. To help us find the right v^, we used the 
epicyclic approximation of Appendix A (see also Maciejew- 
ski & Sparke 1997). In general, consecutive iterates occupy 
the interior of a thick ring, but there is a range of V:^ for 
which the ring generated by a particle starting at that point 
looks almost 1-dimensional. We assume that the real loop is 
nearby, and we need a practical definition of the loop: given 
A'^ (x,v) pairs, we want to know if they lie on a 1-D curve. 
In principle this is not possible: there is a curve through any 
finite set of points. Nevertheless, we are looking for simple 
curves. The ideal choice might be to find the shortest line 
connecting the phase space positions of all the points, and 
to minimize its length as ii^ varies. It is known though, that 
this task (the 'traveling salesman problem') is very expen- 
sive numerically: the number of operations needed to solve 
it grows exponentially with the number of points. Instead, 
we search for the loop by varying in order to minimize 
the width of the ring in the xy plane formed by the iterates 
of the assumed starting point. 

Before minimizing, we adjusted each loop candidate for 
fiattening in the xy plane by finding the points with maxi- 
mum a:-coordinates Xmax and maximum y-coordinates j/max- 
Then we rescaled the y-coordinates, multiplying them by 
J/max/imax, to give a roughly circular ring. This was divided 
into 20 sectors, and the points with the smallest and largest 
radii were found in each sector. The difference of these ex- 
tremal radii gives the width of the ring in each of the 20 sec- 
tors. We define the overall ring width w{y,Vx) as the max- 
imum width, or its average over the sectors. Now, we can 
minimize 'w(y,Vx) with respect to the a;- velocity and get the 
minimum width w{y). Both minimizations proved to be use- 
ful, and the method itself is efficient and fast. In principle, 
we should be able to bring the ring width arbitrarily close 
to zero by increasing the number of iterates and sectors, and 
adjusting the initial velocity Vx- In practice, 20 sectors and 
400 iterates proved to be good and efficient enough. A clear 
minimum is present at the same location for both estimates. 
We assume that the loop is generated by a particle with the 
starting x-velocity for which the width is minimal. The evo- 
lution of one loop thus defined is shown in Figure |^, together 
with the corresponding solution of the epicyclic approxima- 
tion for comparison. The epicyclic approximation remains 
close to the real loop, which justifies using this approxima- 
tion to begin the loop search. For non-symmetric loops, we 
would have to minimize with respect to 2 parameters; for 
example the initial values of Vx and Vy. 
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Figure 2. Evolution of a loop in Model 1 in the frame rotating 
with the main bar, which remains horizontal. The straight solid 
line marks the major axis of the small bar. The epicyclic approx- 
imation is plotted as a solid line, the true loop is dotted. Units 
are in kiloparsec. 



times be large, but the loop family preserves continuity. We 
therefore accepted some fairly shallow minima in the auto- 
matic loop search. Doing so could derail the search from one 
loop family to another, if they come close to one another in 
the phase space (y, Vx)- To clarify when the border between 
different families has been crossed, we searched for orbital 
families twice. Once we did it by following particles start- 
ing on the y-axis, when the bars were both aligned on the 
a;-axis, as described above. We did a second search, starting 
the particle on the j/-axis when the big bar was still on the x 
axis, but the small bar was perpendicular to it. If the search 
algorithm loses the loop, one can go back to the phase plane 
display: identifying the region where the loop was lost often 
helps to answer why it was lost. Examples of this procedure 
will be given in sections 5 and 6. 



Iterates of a given point are not always spread uniformly 
around the ring, but sometimes collect into several distinct 
clumps. This happens when the frequency u with which the 
particle circumnavigates the loop is almost commensurable 
with the potential oscillation frequency uip. If these frequen- 
cies remain close to a lower order resonance, the number of 
clumps is small, and there are few points outside the clumps. 
This may produce an additional spurious minimum in the 
average ring width, but the minimum is narrow and changes 
quickly with the starting y-coordinate when the whole loop 
family is examined (see below), so its character is easily rec- 
ognized. 

In a single bar, some of the periodic orbits are self- 
intersecting; we expect to find self-intersecting loops in dou- 
ble bars as well. Our method of minimizing the ring width 
only in position- or velocity-space obviously breaks down 
there. We could consider the loop in the 4-dimensional phase 
space, where it never intersects itself, at the cost of a con- 
siderably larger computational expense and less clear visual 
examination of the results. Alternatively, we can calculate 
the width for part of the ring only. For example, if the loop 
intersects itself on the a;-axis, these regions can be excluded 
from the calculation of the width. This procedure proved to 
be very fcist and robust, and has been used exclusively. 

The epicyclic approximation provides at most one so- 
lution for the loops and orbits at a given guiding radius. 
Additional loops, not predicted by the epicyclic approxima- 
tion, can be present and may be dynamically important in 
the system of two bars. The only way to find them is to 
search the phase space of the initial conditions {y,Vx) for 
the ring width w{y,Vx) values. In the following sections, we 
call this procedure the phase plane search. Sets of points 
{y,Vx) with low w that group along continuous lines should 
indicate loop families: continuous functions Vx {y) for a range 
of y and Vx ■ The larger the region of low w around the line 
marking a loop family, the more orbits such a loop family 
traps around itself, and thus the more useful 'backbone' for 
building the galaxy it is. 

We searched for loop families by changing the initial 
value y on the j/-axis, and then adjusting the velocity Vx, so 
that the new pair (y, Vx) generates another loop. By analogy 
with periodic orbits in a single bar, we call Vx{y) the char- 
acteristic curve for a given family. In our search we noticed 
that the minimal width of the ring of iterates can some- 



4 MODELS BASED ON OBSERVATIONS 
4.1 The structure of the models 

To minimize the number of free parameters characterizing 
our models, we followed Athanassoula (1992a) in adopting 
analytical forms for the various potential components. Our 
potential consists of a bulge, a disk and two bars. There is no 
halo: in these 2-dimensional calculations its contribution to 
the rotation curve is included in the disk potential. Mass can 
be distributed arbitrarily into spherical and disk-like com- 
ponents, as long as they together mate a realistic rotation 
curve. 

The volume density of the bulge is given by a modified 

Hubble profile 



(4) 



It requires two free parameters: the bulge central density pi, 
and the characteristic radius r^. The surface density of the 
disk is defined by a Kuzmin-Toomre profile (Kuzmin 1956; 
Toomre 1963) 



cT(i?) = ao(l + 7?V-Ro) 



(5) 



like the bulge, this requires two parameters: cto and Ro- Fol- 
lowing Athanassoula, wc adopted the Ferrers ellipsoid to de- 
scribe both bars, with volume density 



Po(l 




if m< 1 
otherwise 



(6) 



where in Cartesian coordinates = /a^ ■\-y^ j}? ■\-z'^ jc? . 
Thus the isodensity surfaces of our bars arc concentric el- 
lipsoids of constant axis ratio, and each bar is described by 
5 free parameters: the central density po, the Ferrers expo- 
nent n, and 3 semi-axes. We adopt prolate bars, and will set 
n = 2, which reduces the number of free parameters to 3. 
A larger index n allows the higher derivatives of the poten- 
tial to be continuous, which in turn reduces the fraction of 
chaotic orbits in the solution. If there are loops supporting 
doubly barred potentials, they are more likely to be found 
for bars with higher n. In Appendix B, we show how the 
potential of a prolate Ferrers bar can be obtained by reduc- 
ing Pfenniger's (1984) solution for a triaxial bar. Both bars 
rotate about the a-axis with respect to the inertial frame; 
their pattern speed, or equivalently the Lagrangian radius, 
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Table 1. 12 parameters for our double-barred galaxy models 

Model 1 Model 2 

disk parameters 

Hmax = V2Ro 20 kpc 20 kpc 

?;niax = 2y'iTGRoao/3°-''^ 164 km s"! 164 km s"! 

bulge and bar parameters 

total mass within 10 kpc A/tot 5 X 10^° Mq 5 X 10^° Mq 
central density p,, + poB + POS 1 x 10^" Mq kpc~^ 4.8 X 10^° Mq kpc~^ 

primary bar semimajor axis ag 7 kpc 6 kpc 

primary bar quadrupole moment Qm 2.25 X 10^" Mq kpc-^ 4.5 X 10^" Mq kpc-^ 

Lagrangian radius of the big bar 1.2ag 1.0a g 

axial ratio of the primary bar ag/bg 2.5 2.5 

axial ratio of the secondary bar ag/bg 2.5 2.0 

ratio of bar semimajor axes ag /as 0.6 0.2 

ratio of bar masses Mg/Mg 0.6 0.15 

pattern speed of the secondary bar Qg 42 km s~^kpc~^ 110 km s~^kpc~^ 



constitutes a 4th free parameter for each bar in our case. 
Thus the whole system is modeled with 12 parameters. 

The parameters for the two models considered in this 
paper are given in Table 1. The potentials of the bulge, disk, 
and the big bar are chosen on the basis of Athanassoula's 
(1992a) 'standard model', constructed by Wozniak (1991) 
from a set of potentials inferred from observations of 16 SBO 
galaxies and one SBa. The disk parameters are the same as 
in Athanassoula's model. The parameters of the secondary 
bars were chosen on the basis of the observations listed in 
the next subsection. In Figure ^, we show the rotation curve, 
calculated from forces along major and minor axis of the bar 
at bars aligned, and the curve of f2 — k/2 in the azimuthally 
averaged potential for both models. 

Model 1 was designed to explore how far the secondary 
bar can extend in a doubly barred potential. It is the same 
model as the one used by Maciejewski & Sparke (1997). The 
parameters of the axisymmetric part match Athanassoula's 
average observed values almost exactly. But the main bar is 
weaker than the observed average, although falling within 
the range found by Wozniak (1991). 

In Model 2, the quadruple moment of the primary bar, 
its Lagrangian radius, and its axial ratio are the same as 
in Athanassoula's model. Since we use n=2 Ferrers bars, 
whose density drops faster with radius than in the n=l bar 
modeled by Athanassoula, we extended the main bar all the 
way to its Lagrangian radius at 6 kpc, so the mass distri- 
butions are roughly similar. We increased the total central 
density to a value twice than that in Athanassoula's model, 
so the secondary bar comes close to having an ILR (Fig. |^. 
Of the 16 galaxies in the Athanassoula & Wozniak sample 
(Wozniak 1991), only two have such a big maximum of the 
Q — k/2 curve, but the other parameters for the axisymmet- 
ric part fit these observations well. Both bars are stronger 
here than in Model 1, though still weaker than Wozniak's 
observed average. 

4.2 Observed parameters of doubly barred 
galaxies 

Friedli (1996) reviews available observations of possible dou- 
ble bars and estimates their basic parameters. We follow his 
notation here, in which three parameters are defined: 



a) Model 1 



b) Model a 
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Figure 3. The rotation curve R d^/dR) and its components 
along the major axis (top panels) and minor axis (central panels) 
for a) Model 1 and b) Model 2. The bottom panels show the 
angular frequency curves in the azimuthally averaged potential 
that determine the positions of resonances in the big bar (B) and 
small bar (S). We marked the inner and outer ILRs (ilLR and 
oILR respectively), the corotation (CR), and the Outer Lindblad 
Resonance (OLR). 



(i) the maximum ellipticity in the barred region e™^" — 
maxa |1 — fe/a|, where a and b are the semimajor and semimi- 
nor axes of the fitted ellipses (note that this is the maximum 
ellipticity of the total hght distribution, not of the bar only), 

(ii) the length ratio of the bars /3 — Ig/lg, where lg,lg are 
the semi-major axes of the isophote ellipse fits at maximum 
or minimum ellipticity that corresponds to the main and 
secondary bar, respectively, 

(iii) the luminosity ratio of the bars 7 = Lg/Lg: the ratio 
of total luminosities within the isophote of maximum (min- 
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Table 2. Observed parameters of doubly barred galaxies (Friedli 
1996 / Mulchaey et al. 1997) 





average 


min 


max 




7.2/4.9 


3.7/1.83 


18.0/8.8 


7 


4.0/- 


2.0/- 


7.5/- 




0.34/0.30 


-/0.2 


-/0.43 




0.55/0.41 


-/0.2 


-/0.6 



imum) ellipticity (note that thus defined, the inner bar's lu- 
minosity is strongly influenced by the light from the bulge). 

Wozniak et al. (1995) discuss how the last two parame- 
ters above depend on whether they are derived at minimum 
or maximum ellipticity. Mulchaey, Regan & Kundu (1997) 
present another, smaller but rather unbiased sample of dou- 
bly barred galaxies seen in K-band images. Average values 
at minimum ellipticity of parameters defined above for both 
samples (13 and 6 doubly barred galaxies respectively) are 
given in Table 2. There is an observational bias towards 
lower values of /3 and 7, because such systems are most eas- 
ily found. Also e™'"' for the small bar may be underestimated 
because of finite spatial resolution. 

4.3 Model parameters versus observational 
constraints 

To make a crude comparison between these models and the 
observed double bars, Erwin (1998, private communication) 
performed ellipse fitting on both our models in face-on po- 
sition for a relative bar orientation of 60°, and derived the 
values of /3, 7, e™'"' and eg"'^. Instead of surface bright- 
ness, the surface density was used. Since the contribution 
of the dark halo in our models is included in the disk mass 
distribution, the halo appears luminous here. It is not im- 
portant though, since we are interested in the central parts 
of the galaxy, where the luminous matter is dominant. We 
calculate the luminosity of each bar as an integral over the 
interior of the ellipse with maximum or minimum ellipticity. 
We consider the parameters derived at the maximum ellip- 
ticity (max-e) to be better defined, because in real galaxies 
the radius of minimum ellipticity (min-e) at the main bar 
depends on the inclination angle. Table 3 gives the derived 
parameter values. 

The /3 and 7 parameters for Model 1 derived with both 
'min-e' and 'max-e' methods are too small compared to the 
observed range (Table 2). Thus the secondary bar in Model 
1 is too large; despite that, we explore this model to check 
if such a configuration is theoretically possible. 

In Model 2, we see well defined ellipticity maxima, with 
a minimum between them. Note that both 'max-e' and 'min- 
e' methods give bar sizes smaller than the major axis of 
the Ferrers bar: in Model 2 the primary bar's half-length 
is estimated to be 3.58 and 5.58 kpc respectively, when the 
semi-major axis of the Ferrers primary bar is 6.0 kpc. The 
measured values of the half-length of the secondary bar are 
0.62 and 0.92 kpc respectively, when we set as=1.2 kpc. The 
/3 and 7 parameters in 'min-e' and 'max-e' estimations have 
very similar values. Note that the mass in the inner parts 
of the galaxy is dominated by the bulge, and this makes 7 



Table 3. Model parameters derived from isophote fitting 





Model 1 


Model 2 




min-e max-e 


min-e max-e 




3.07 1.98 


6.08 5.79 


7 


2.38 1.65 


3.55 3.58 


^max 

'^s 


0.118 


0.159 




0.182 


0.428 



much smaller than the mass ratio of the bars. For Model 2, 
the f3 and 7 values derived by the ellipse-fitting method cor- 
respond very well to the parameters in the observed samples 
(Table 2). 

The only parameter of Model 2 that does not conform 
to the observed range, is the 'measured' maximum ellipticity 
of the secondary bar. This is only half as large as observed 
average, because the central parts of the model are domi- 
nated by our strong spherical bulge. It is very likely then, 
that distinct secondary bars are best observed against weak 
bulges, which cannot induce an ILR in them, unless the bars 
themselves are highly centrally concentrated. 

From Figure H we see that the corotation of the sec- 
ondary bar in Model 2 (as=1.2 kpc) is at 2.3 kpc, well out- 
side the bar. Such a setup was forced after Model 1 failed to 
provide orbits supporting outer parts of the secondary bar 
— see the following section. The slowly rotating secondary 
bar seems to be consistent with the most successful numeri- 
cal model by Friedli & Martinet (1993, Set II). They observe 
a secondary bar of 2.5 kpc size with corotation at 3.5 kpc, 
which makes r_t/a ~ 1.4, compared to 1.9 in our case. Thus 
the slowly rotating secondary bar appears to be justified by 
numerical models. 



5 EXAMPLES OF LOOP SEARCH METHODS 
5.1 Loop families in a single bar 

To illustrate the methods which allow us to find loops in 
doubly barred galaxies, we show how the loop search is per- 
formed for a single bar. Here, we can choose arbitrarily the 
period of integration (we chose it to be very small) used to 
generate consecutive iterates of a given point. If this point 
lies on a periodic orbit, all its iterates will also he on this 
orbit: in a single bar our loop finding technique locates the 
closed periodic orbits. If the orbit is stable, nearby particles 
should generate points populating a ring around this orbit. 

We performed the search for Model 2, with the mass 
of the inner bar redistributed into the axisymmetric bulge, 
so that the central density and the total mass defined in 
Table 1 remain unchanged. Figure ^ illustrates the phase 
plane search defined in section 3.2. Two dark arches sur- 
round the main stable orbital families; in the notation of 
Contopoulos & Papayannopoulos (1980), these are the xi 
family, elongated parallel to the bar, and the X2 family, ori- 
ented perpendicular to the bar. Thick dark arches, which 
overwhelmingly dominate the phase plane, mean that many 
quasi-periodic orbits are trapped around these periodic or- 
bits. The unstable X3 family is absent from this plot, because 
particles are not trapped around it. 
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Figure 4. a) Phase plane search for particles in Model 2 without 
the secondary bar, which start on the minor axis of the bar. The 
colour reflects how strongly the first 400 iterates of a particle with 
given starting values {y,Vx) are confined to a 1-D curve. Darkest 
areas mark the near-loop solutions (smallest ring widths to), b) 
The characteristic curves for Model 2 without the secondary bar. 
The diagrams of minimal ring widths w(y) along each family are 
attached. 



The characteristic curves Vx(y) for the xi and X2 fam- 
ilies are shown in Figure together with the minimum 
values of the ring width w defined in section 3.2, describ- 
ing how close a set of iterates of the starting point given 
by {y,Vx) falls to a 1-D curve. The smooth appearance of 
the widths assures us that we follow a single loop family. In 
particular, when the orbital family X2 disappears at large y, 
the ring widths increase rapidly. In §6, we will compare these 
orbits to the loops in the doubly barred Model 2, which is 
close to being self-consistent. 

5.2 Loop families in a generic doubly barred 
potential: Model 1 

Here we discuss the loop families in Model 1, the first of 
our doubly barred potentials. It is more general than Model 
2, because the bars are not in mutual resonance, and the 
secondary bar is relatively large. Both bars are rapidly ro- 
tating: they extend to ~ 85% of their Lagrangian radii. To 
enable a comparison of the calculated loops to those from 
an epicyclic approximation, we made the bars weaker than 
in Athanassoula's (1992a) standard model. Both bars have 
two ILRs: those of the main bar are at 0.32 kpc and 2.9 kpc, 
and of the small bar at 0.6 kpc and 1.3 kpc. 

Figure ^a presents a phase plane search for loops which 
are symmetric with respect to both x- and j/-axes at the 
time when the bars are aligned on the x-axis. As in Fig- 
ure ^a, two dark arches surround the stable loop families, 
though the top arch is broken in the middle. Extending the 
notation of orbital families, we call the loop family that cor- 
responds to the X2 orbital family of the main bar, the X2 
loop family. It appears clearly in Figure along with a 
wide dark stripe of near-loop solutions, which shows that 
particles can easily be trapped around this family. The xi 
loop family corresponds to the x\ orbits in the outer bar. 
We see the inner part of the xi loop family clearly out to 
the radius where the bifurcation giving rise to the X2 fam- 
ily takes place. The family abruptly discontinues above this 
point [vx «250 — when we refer to positions in Vx(y) dia- 
grams in this and the following sections, the length units are 
kpc, and the velocity units are km s~^). The inner part of the 
xi loop family remains inside the inner ILR of the main bar. 
The xi loop family can be detected again for 0.4 < y < 2.0, 
between the two ILRs of the main bar. But here it is not 
accompanied by a wide dark region of near-loop solutions, 
indicating that it may be not able to trap many particles 
around itself. Our plane search finds no clear xi loop fam- 
ily in the region 2 < y < 3, close to the outer ILR of the 
main bar; there may be no loops there, or the resolution of 
our search may be too low — it shows only a grey region of 
near-loop solutions in this y range. For larger y this family is 
seen again, together with the accompanying dark grey band 
of near-loop solutions. The discontinuity of the xi family at 
large Vx seems to be common in our double-bar potentials; 
therefore we divide this family into an inner part xu inside 
the Vx{y) maximum and an outer part xio outside it. 

In addition to the xi and X2 loop families, there is an in- 
teresting dark line marked in Figure ^ by A, which departs 
from the inner xi family in the lower left corner of the phase 
plane {y « 0.13, Vx ~ 60). It is examined separately below. 
No other loop families can be clearly seen in Figure In 
particular, there seem to be no stable loops corresponding 
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Figure 5. a) As Fig.|4|a, but for particles in Model 1, which start 
on the minor axis of the bar when the bars are aligned, b) The 
characteristic curves for Model 1. Loop families and resonances 
are labeled: xi; and xio are the inner and outer parts of the xi 
family, as described in the text. The thickness of the line corre- 
sponds to the ring width: the loops with clear minima in width 
are represented by the thickest lines, the dotted lines mark areas 
where loops are no longer well defined. The labels A and B refer 
to loop families described in section 5.2.3 



to X2 orbits of the small bar, although this bar has a pretty 
wide zone between inner and outer ILRs. 

All the loop families that we have found for Model 1 are 
displayed in Figure they are described in detail in the 
rest of this section. The thick lines correspond to well defined 
loops of small ring width; thin lines represent loops that are 
still well defined as single minima, but with considerable 
ring widths. The broken lines correspond to regions where 
loops were not well defined, mostly because at fixed y the 
ring width had multiple minima with Vx, or the minimum 
value of w was rapidly growing with y. 



5.2.1 The X2 loop family 

We were able to follow the X2 family as a continuous min- 
imum of the ring width from (y=G.A, Wa; = 130) to (j/=2.2, 
Wa; = 140), which remains inside the zone between inner and 
outer ILRs of the main bar. The X2 loops still can be found 
in the range 2.2< y <2.8, but the X2 family is broken in this 
region. Our attempts to follow it fail, because multiple ring 
width minima compete with one another as y changes. 

A representative set of the X2 loops in the frame rotating 
with the primary bar is shown in Figure ^a. Here we confirm 
the hypothesis that we developed on the basis of the epicyclic 
approximation (Maciejewski & Sparke 1997): the outermost 
X2 loops correspond to the X2 orbits of the main bar and 
remain roughly perpendicular to it, while the inner loops of 
this family follow the secondary bar in their motion, and 
correspond to the xi orbits of the secondary bar. However, 
we cannot trace the X2 loops all the way to the center, while 
the xi orbits of a single bar generally do extend to the center. 

Throughout this paper we assume that our 2 bars are 
rigid and uniformly rotating. But the structure of the X2 
loops in Figure ^ tells us that it is a crude approximation 
only: the secondary bar, as outlined by its orbital (loop) 
structure, pulsates and accelerates as it moves through the 
main bar. The loops supporting the small bar are much 
thicker when it is parallel to the main bar: the minor axes of 
the inner loops are up to 1.8 times as large as when the bars 
are perpendicular. So if a secondary bar is self-consistently 
made out of particles trapped around these loops, then it 
should be thinner while perpendicular to the main bar. Also, 
the X2 loops lead the small bar as it moves from the parallel 
to perpendicular position, and trail behind it when it moves 
back to align with the main bar. One can conclude that in a 
self-consistent model, the small bar should decelerate when 
going out of alignment with the main bar, and accelerate 
when on its way back to the aligned position. These con- 
clusions will be confirmed by the loop structure of Model 
2. 

In Model 1, the X2 loop family extends only to x=2 
when the bars align on the x-axis, and only the loops ex- 
tending roughly to x=1.5 there follow the small bar in its 
motion. But the secondary bar itself extends to as=4.2 kpc, 
so there are no loops which support the outer part of the 
secondary bar. By analogy with argument of Contopoulos & 
Papayannopoulos (1980) for the xi orbits in a single bar, we 
conclude that a self-consistent secondary bar cannot extend 
outwards beyond the region where the X2 loops follow its 
motion relative to the large bar. The X2 loops populate the 
region between the ILRs of the large bar. So we conclude 
that the major axis of the small bar is unlikely to be larger 
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Figure 6. a) The X2 loops in Model 1 in the rotating frame of the 
large bar, which remains horizontal. The straight lines at opposite 
sides or corners indicate the position angle of the small bar; time 
increases down the page, b) The xi loops at the same relative bar 
phases. The inner- and outermost X2 loops are marked by dotted 
lines; three xii loops are displayed inside them. Units are in kpc. 



than the outer extent of the X2 orbital family of its main 
bar host along its major axis, unless the strength of the sec- 
ondary bar relative to the main one is large enough to cause 
major changes in the loop structure. 



5.2.2 The xi loops 

The inner xi loops (the xu subfamily) show a very regular 
behavior in Figure Ha, where they are seen as a thick line in 
the lower left corner. The Vx{y) curve is continuous through- 
out the whole region where these loops were detected: from 
the very center of the potential to y— 0.223. At this y, the 
a;-velocity increases rapidly with y and the loop family is 
lost in Figure The xu loops are almost circular and are 
represented by three innermost curves in Figure 

Unlike the xu loops, loops belonging to the xio sub- 
family (outer xi loops) do not define a continuous function 
Vx(y) in the plot of Figure The loop searching code re- 
peatedly lost track of the minimum after covering only a 
short y range. For some y ranges, two comparably deep and 
wide minima were present, indicating that two loops exist 
at the same y but with different starting velocities Vx. The 
appearance of the Xio loops is shown in Figure |^. It is not 
regular, and these loops cover a very small part of the pri- 
mary bar. Therefore it appears unlikely that the primary 
bar can be built in a self-consistent manner in Model 1. 



a 
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Figure 7. a) Evolution of the representatives of the loop fam- 
ily marked by A in Figure ^ displayed in the rotating frame of 
the large bar, which remains horizontal. The straight line indi- 
cates the position angle of the small bar. Only loops for the bars 
coming out of alignment are displayed: as the bars come back to 
alignment, the loops are mirror reflections of these, b) Represen- 
tatives of the family marked by B in Figure displayed in the 
same way as those in a). Labels are in kpc. 



5.2.3 Other loops 

Although the secondary bar has an ILR, we found no loops 
corresponding to its X2 orbits. Instead, we found a family 
of loops indicated by the set of dark pixels in Figure Bp, 
(marked there by A), that departs to the right from the 
thick band marking the xn loops. Although this family is 
not surrounded by a dark band of particles trapped around 
these loops, the minimum ring width is well defined, and 
loops representative of this family are displayed in Figure 
. When the two bars are aligned on the a;-axis, the loops 
are perpendicular to them. But the loops precess in a sense 
opposite to the small bar, aligning with it when it makes an 
angle 7r/4 with the main bar, and becoming perpendicular 
to the small bar again when it in turn is perpendicular to 
the main bar. Thus the figure of the loop moves retrograde, 
but loop particles remain on prograde orbits with respect to 
either bar. The inner loops become rounder, and this family 
joins the xu loops at the point y—0.13, Vx~57 (Fig. |^). 
This family obviously does not correspond to the X2 orbital 
family in the small bar, and it does not support either bar. 
Nevertheless, it is consistently present in further modifica- 
tions of Model 1, and in Model 2 described below. 
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Just as we can have multiply periodic orbits in a sin- 
gle bar, we can have higher order loops in double bars. If 
particles come back to the same curve only every other bar 
alignment, then after an odd number of alignments they 
lie along one curve, and after an even number, on another 
curve. Thus a single particle generates two closed curves in 
this case. One can easily extend this to curves generated 
by a single particle returning to the same curve after every 
A'^ alignments of the bars. We found many such higher or- 
der loop families in our search, one of which is presented in 
Figure both curves constituting each loop there precess 
retrograde, with a period twice the the relative period of 
the two bars. In Figure ^b, this family (marked by B) runs 
parallel to the abovementioned family, when the larger out 
of two y-axis crossing points is considered to be a loop gen- 
erator. Obviously, this family does not support any bar, and 
it is very unlikely that any of the higher order loop families 
can do it. Although we expect that abundant and complex 
higher order loops are present, we did not search for them. 



6 APPROACH TO A SELF-CONSISTENT 
DOUBLY BARRED MODEL 

Louis & Gerhard (1988) have shown, and the previous sec- 
tion illustrates, that bars which make up a self-consistent 
double-bar system must distort slightly as they rotate 
through each other. Thus our assumed potential, with two 
rigidly rotating bars, cannot be fully self-consistent. Here 
we are looking only for a model that may be close to being 
self-consistent, in which stable regular orbits give support 
to both bars in their motion. A real bar would have some 
chaotic orbits as well, which we do not consider here. 



6.1 Choice of parameters leading to a nearly 
self-consistent model 

A system similar to Model 1 could not be built in a self- 
consistent manner: its X2 loops do not extend to the end of 
the secondary bar, and the xi loops cover only a small frac- 
tion of the primary bar. Making the secondary bar stronger 
may help to extend the X2 family, but it would destroy the 
xi loops, which are already quite disrupted by the small bar 
in Model 1. A different solution has to be found. 

First, we modified Model 1 by increasing the a/b axial 
ratio of both the bars from 2.5 to 4.0, while the product ab 
of the axis lengths remained constant. The xi loop family 
virtually disappeared, except for the inner xu loops. The 
X2 family broke into two well separated parts, and it may 
be difficult to support the secondary bar as well. On the 
basis of our limited search, weakness of the xi family seems 
to be a general feature of models with a relatively large, 
fast-rotating inner bar, like Model 1. 

The new model, hereafter Model 2, is meant to be as 
close to Athanassoula's (1992a) standard model as possible. 
The quadruple moment of the primary bar is twice as large 
as in Model 1. From Model 1, we know that the secondary 
bar should not be too big. First we tried a long, but not mas- 
sive {as = 0.4as, Ms /Mb = 0.2), secondary bar. The phase 
plane search showed that the x\ family, which should sup- 
port the longer bar, is practically absent for y >0.9, while the 
bar length is as =6. The X2 family is broken into two parts. 



both surrounded by regions of trapped semi-periodic orbits. 
This model may have a nearly self-consistent secondary bar, 
but cannot support the primary bar, and therefore has to 
be rejected. 

Then we tried a smaller, but more massive secondary 
bar (as = 0.2as, Ms/Mb = 0.3). The loops from the xi 
family extend far enough in radius (1.1 < y < 4) to support 
part of the primary bar. The inner part of the xi family is 
quite extended, though: it suppresses the X2 family, which 
in turn appears irregular and disrupted on the phase plane. 
In this model, part of the primary bar can be supported by 
loops, but the secondary bar is likely to be far from self- 
consistent. 

The two attempts above suggest that our secondary bar 
is still too massive, or too large. Observations presented in 
§4 confirm this suspicion: the ratio of semimajor axes of 
the bars in the model, is beyond the observed range for our 
model with as ~ 0.4as, and at the lower end of the ob- 
served range for the model with as = 0.2as. Guided by the 
values listed in Table 2, we decided to make the secondary 
bar smaller and less massive. The ellipticity of the observed 
secondary bars is considerably smaller than that of the pri- 
mary ones; therefore we set the axial ratio of the secondary 
bar at bs/as = 0.5. The small bar is only 20% as long as the 
primary bar and has only 15% of the big bar's mass. Since in 
Model 1 the loops supporting the small bar failed to extend 
all the way to its corotation, in Model 2 we set the pattern 
speed so that the Lagrangian radius of the secondary bar 
is at 1.9as. This also puts corotation of the small bar at 
the outer ILR of the big one, as suggested by Tagger et al. 
(1987). 



6.2 Loop families supporting a nearly 
self-consistent Model 2 

We constructed Model 2 with a particular interest in loops 
supporting the potential. We want to answer the question, 
is it possible to construct a model close to this one in a 
self-consistent way? Our first task is to survey the structure 
of x\ and X2 loop families, which we expect to support the 
primary and secondary bar, respectively. The results of the 
phase plane search for Model 2 are displayed in Figure |^. 
Both x\ and X2 loop families are present there, with trapped 
orbits around them indicated by the dark areas. Note how 
similar this figure is to Figure Ha, in which a phase plane 
search for Model 2 with the main bar only was presented. 
Unlike in Model 1, this plane search suggests that both bars 
are likely to be supported by stable loops. Figure ^ shows 
the characteristic curves and the variation of minimum ring 
width w(y) along each loop family. We find an additional 
loop family, marked there by br , supporting part of the main 
bar. 



6.2.1 The X2 family — loops supporting the small bar 

We were able to find loops belonging to the X2 family for 
0.17 < y < 2.39. Just before the X2 family is lost at small y, 
the velocity drops rapidly. At the high-y end of the X2 family, 
the ring widths increase making it harder to find the loop, 
before the loops disappear (Fig.Hb). This family has two in- 
teresting features. First, around j/=0.406, the characteristic 
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Figure 8. a) Phase plane search for Model 2 with two bars ini- 
tially on the i-axis. b) The characteristic curves for Model 2 at 
bars aligned. The diagrams of minimal ring widths 'w{y) along 
each family are attached. For clarity, some characteristic curves 
and the corresponding widths are marked in dotted. Note that 
the widths change much more abruptly with y than in Figure M. 



curve is discontinuous, and there is a small y range where no 
loops have been found. This corresponds to the white stripe 
crossing the X2 arch in the phase plane search of Fig. The 
minimal ring width w{y) increases steadily as this disconti- 
nuity is approached from both lower and higher y's, before 
the width experiences a sudden jump. A detailed search for 
loops in this y range found no simple smooth invariant 1-D 
curve; but for a wide range of Vx, the particles are confined 
to well-defined rings. No higher order loops were detected 
in this region. It is possible that a similar small discontinu- 
ity occurs at ^=1.86. Here, the ring width jumps an order 
of magnitude above its normal level over a narrow j/-range, 
where we find no simple smooth 1-D curves. Again, for a 
wide Vx range we observe that the particles are confined to 
rings there. 

At 1/ ~ 1.08, the Vx{y) curve for the X2 loops splits 
in two; over a small range, there are two loops for each y 
value. The accompanying plot of ring widths shows that 
when approaching this region from higher y, the width of 
loops increases rapidly, and the searching code effectively 
loses the loop family. The minimum followed by our contin- 
uous search is no longer the global minimum, and another 
minimum is found which generates the continuation of this 
family to smaller y. 

The overall appearance of the X2 loops in Model 2 is 
presented in the central panels of Figure ^. The X2 loops 
with y ^ 0.5 rotate with the small bar and extend out to 
about 0.85-0.9 of its semimajor axis. There the density in 
an n=2 Ferrers bar has dropped to 4% of its central value. 
Thus we conclude that the secondary bar is well supported 
by X2 loops: even if the family breaks at j/ = 0.4, there are 
'near-loop' solutions, and particles trapped around these can 
serve as the building blocks of the inner bar. Here we can 
also confirm our conclusion drawn on the basis of Model 1, 
that the secondary bar pulsates and accelerates, although 
these effects are smaller in Model 2. 



6.2.2 Loops supporting the larger bar 
families 



the x\ and bj 



Figure |8p shows two disconnected parts of the xi family, 
which were detected with the search method described in 
section 3.2: the xu loops (inner xi) a.t y < 0.09 and the xio 
loops (outer xi) at j/ > 0.5. 

The xii loops show a rapid growth in the a;-velocity with 
y, and for Vx >120, this loop family continues for decreasing 
y with velocities still growing. Therefore the xu loops are 
double- valued in this range. They cannot be traced beyond 
the point Vx=310, y=0.024, corresponding to a dark pixel in 
Figure ^a. We explore the question of connection between 
the inner and outer xi families later in this section. 

The xio loops are not as well defined as the xu and X2 
loops — the minimum ring widths are larger. Near j/=0.65 
and 0.85, they rise to unacceptable values of 50% and 30% of 
the ring size, but no discontinuity in the loop family tracing 
is seen. We will return to this feature later in this section. 
In general, for the xio loops there is no clear single w{vx) 
minimum for a given y, but rather several less deep minima, 
of which one usually dominates. This is similar to what the 
X2 loops experience near y = 1.08. For xio loops, two minima 
of roughly equal depth persist for the y range from 1.85 to 
2.4, and are displayed as two parallel lines in Figure pp. In 
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Figure 9. An overview of loops in a nearly self-consistent Model 2 of a doubly barred galaxy. The loops are displayed in the reference 
frame of the main bar at different relative positions of the bars. The outer bar is outlined by the dashed line and the major axis of 
the secondary bar is marked by a straight line on outer panels. The sequence follows along outer panels clockwise, with central regions 
magnified on inner panels, where the secondary bar is outlined by the dashed line. Units on the axes are in kiloparsec. 



this region, the ring widths increase with y for the loops 
marked with a solid line, and this branch of the xio loops 
is eventually lost at y=2.4. At that y, the widths on the 
other branch, marked with a dashed line, become small and 
behave in a very regular way. This branch is lost in turn 
at j/=2.65. At 2.57< y <3.22, another branch of the xio 
loops exists; this is the branch extending to the highest y in 
Model 2. We consider the solid line extending from y=0.5 
to 2.4 to be the main part of the xio family, and we call the 
two abovementioned branches the xio family additions. 

The dotted line parallel to the xio family and persisting 
for the y range from 0.54 to 1.4 (Fig. |^), marks another 
loop family found in Model 2, which we name 6t for their 
'bow-tie' appearance. Below 1 kpc, the ring widths for the 
6t loops are smaller than those of the xio loops, and the 
widths behave in a more regular way. 

Figure ^ displays representatives of all loop families 
found for Model 2, as they change while the two bars ro- 
tate relative to each other. The two outermost loops belong 
to the two outer additions to the xio family. The next four 



loops inside them belong to the main part of the xio fam- 
ily, and they undergo interesting oscillations when the two 
bars rotate through one another. The loop that is inner- 
most when the bars are aligned (top left panel) develops 
cusps, and becomes much rounder: it ends up as a second 
outermost loop when the bars are orthogonal (bottom right 
panel). This leaves us with no a^io loops supporting the inner 
part of the main bar when the bars are orthogonal (bottom 
right panel). But there is one more loop inside the above 
mentioned set of four, which remains parallel to the main 
bar. It intersects itself at bars aligned, and belongs to the 
6t family marked in Figure]^ by a dotted line. Loops from 
this family take a boxy shape at bars orthogonal and fill the 
central part of the main bar, thus supporting it in addition 
to the xi loops. 

The xii loops are displayed against the X2 loops in the 
4 inner panels of Figure y. The double- valued Vx{y) curve 
for these loops means that when the two bars are aligned, 
two loops from this family cross the y-axis at the same point. 
One set of these loops is close to circular and can be detected 
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all the way to the potential center, the other one is elongated 
and gets longer with decreasing y, up to a semimajor axis 
of 0.5 kpc for i/=0.024, where the family disappears. These 
narrow loops get stretched when the bars are orthogonal, 
but their semimajor axes roughly do not change, and they 
remain well aligned with the primary bar. Note that the 
xii loops are narrowest at bars aligned, unlike the X2 loops 
which are narrower at bars orthogonal. Figure ^ also shows 
that in a potential of two independently rotating bars, sets of 
loops supporting the bars can penetrate through one another 
while each follows a different bar. If there are enough stars 
trapped on the xu loops, such a galaxy may appear as having 
a triple bar, with innermost and outermost bars aligned. 

6.2.3 Loop search with bars orthogonal 

One can get a more complete insight into the loop fami- 
lies in double bars by performing an additional search, with 
particles starting on the minor axis of the main bar (the 
j/-axis), moving perpendicular to that axis, when the bars 
were orthogonal. Phase plane search of Figure indicates 
both the x\ loops (the upper arch) and X2 loops (inner arch 
and plateau). A close-up of the upper arch is shown in Fig- 
ure hOp: many bright strips running from the top-left to the 
bottom-right indicate areas where the iterates are no longer 
confined to a ring. The xi family is broken when they cross 
its dark arch. The structure at such intersections is simi- 
lar to the features common in characteristic diagrams for 
orbits in a fixed potential, usually called the 4/1-like gaps 
{e.g. Contopoulos & Grosb0l 1989). It can be best seen at 
y = 0.2, Vx = 460 in Figure px|b. Around y = 0.6, there is a 
white strip which separates two dark strips: one surrounding 
the hr family and originating at the top-left corner, and the 
other surrounding the x\o loops originating on the right. It 
looks like this white strip overlaps with what would be an 
xi orbital family in a single bar there (see Fig. ^), and it 
is unlikely that stable loops can be found there. 

Comparing the loop families found with particles start- 
ing at the bars orthogonal with those found by the origi- 
nal search method earlier in this section allowed for a more 
complete picture of loops supporting the main bar. Figure 
[ll] shows the characteristic diagrams for these loops, both 
at bars aligned and orthogonal. The xio loops, which occu- 
pied the region from 0.22 to 2.4 kpc in y at bars aligned, 
are compressed into a 1.4 < y < 2.8 range at bars orthog- 
onal. In this configuration, we find no xio loops at j/ < 1.4 
to support the main bar. When the bars are orthogonal, at 
1.4< y <2.2 this family develops an interesting zig-zag on 
the characteristic curve (Fig.|ll|b). It means that the inner- 
most loops at bars aligned get pushed away at bars orthogo- 
nal and end up as the outermost ones. We gave an example 
of such a loop in section 6.2.2. On the diagonal of the 'Z' 
the trend reverses, and loops pulsate without changing their 
order. Eventually, the trend reverses again (the top of the 
'Z'); the loop which is outermost when the bars are aligned 
sinks inside, and is innermost at bars orthogonal. At the 
points {y =1.4, w^=308) and {y—l.&, =250), where the 
xio characteristic curve reverses, there is a set of rings that 
for a short range continues in the direction before reversing, 
but the ring widths increase rapidly. 

In some regions of the lower branch of the 'Z' at bars or- 
thogonal, loops do not transform to points on the xio charac- 
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Figure 10. a)As Fig.g but for a particle starting on the minor 
axis of the main bar when the bars are orthogonal, b) Close-up 
of the top of the 'xi arch' from a). 
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Figure 11. Two diagrams to compare the characteristic curves of 
loop families supporting the main bar in Model 2 at bars aligned 
(top) and orthogonal (bottom). 



teristic curve at bars aligned monotonically, but get spread 
over a large y range. The transformed curve is disrupted 
at several places, which correspond to high ring widths for 
loops found at bars aligned (Fig.^). Some real gaps in the 
stable loop families likely occur there, which are not just 
artifacts of our method. Thus the 'high width' found with 
our method may suggest that the loop sequence is broken 
and we may encounter part of a different family with a short 
stability range. 

The br family and the two disconnected Xio branches 
at high y in Figure [ll| a reverse order at bars orthogonal 
(the innermost loop becomes the outermost one). The two 
outer additions to the xi family may be connected by bridges 
similar to the diagonal of the 'Z' in the main xio subfamily: 
at y = 2.6, = 170, the first addition to the xio loops 
changes order and aims towards the second addition, but 
we were not able to find a continuation for this set of loops. 

We tried to find loops that could connect the xu and 
xio subfamilies. Our search at bars orthogonal resulted in 
finding loops represented by an arch at y—0.2, =450 in 
Figure |ll]b. It may be a continuation of the xu subfamily 
extending to the galaxy center, but no connection of those 
two has been found. Several gaps in the 'xi arch' there seen 
in the phase plane search (Fig.icjb), suggest discontinuities 
in the xi family. The xn arch in Figure |ll| a most likely has 
no connection with xio loops either, since the 6t loop family 
seems to separate them. 



6.2.4 Summary of loops supporting the main bar in Model 
2 

The main xi family in Model 2 consists of 2 subfamilies with 
some additions: 



(i) the inner xn subfamily originating at (y = 0, fx = 0), 

(ii) the addition to the inner xu loops marked by the arch 
in Figure 

(iii) the main xio subfamily: there is no continuous tran- 
sition from this subfamily to the xu subfamily, 

(iv) two outer additions to the xio subfamily. 

There is an additional family, the 6t or bow-tie loops, but 
it still leaves a gap without a continuous set of simple stable 
loops for 0.5< y <1.4 at bars orthogonal. There are some 
small loop subfamilies found at 0.7< y <1.5 at bars orthog- 
onal, but they are discontinuous. There are absolutely no 
loops at 0.5< y <0.7 — there is a white strip of diverging 
iterates overlapping that part of the dark 'a;i arch' on the 
phase space diagram in Fig.^|b. 

In conclusion, we would say that the xi and br loops 
support most of the main bar at bars aligned. The xio loops 
are strongly infiuenced by the motion of the secondary bar. 
At bars orthogonal, they get rounder and group in the outer 
parts of the main bar. We were not able to find any xi loops 
supporting the inner part of the main bar when the sec- 
ondary bar is perpendicular to it - the br loops fill up the 
center there, but a small gap between these two families re- 
mains (FigJol). In section 6.2.1 we showed that the secondary 
bar is fully supported by the X2 loops. 



7 DISCUSSION, CONSEQUENCES FOR GAS 
FLOWS 

The main goal of this work was to explore whether a galac- 
tic stellar disk with two non-equal bars rotating at different 
pattern speeds is dynamically possible. The 'backbone' for 
the orbital structure of this time-periodic system is provided 
by stable loops. These are closed curves, such that particles 
placed on them and moving in the potential of this sys- 
tem will return to the original curve after two bars have 
come back to the same relative orientation. A nearly self- 
consistent double bar should have stable loops elongated 
with each bar over appropriate radial range. In section 6, 
we gave an example of a doubly barred system where stable 
loops existed with the right shape, over the right range in 
radius, to support both the inner and the outer bars. 

To create conditions most favorable to the existence of 
stable loops, we selected a smooth potential by choosing a 
high Ferrers index of the bars {n — 2). We also minimized 
the number of resonances by putting the corotation of the 
small bar at the outer ILR of the main bar. This overlap- 
ping of resonances, although not essential for the existence 
of self-consistent double bars, may be natural in dynamical 
systems (see Tagger et al. 1987). Our secondary bar was slow 
rotating, with corotation well beyond end of bar. To have 
the same resonant condition for a fast-rotating secondary 
bar, the inner bar would have to extend to the outer ILR 
of the main bar. In such a case, all the X2 loops should fol- 
low the secondary bar in its motion, a configuration that we 
were not able to create. On the basis of our models, we be- 
lieve that for double bars we may not be able to have both 
resonance overlapping and a secondary bar extending to its 
corotation. 

The morphology of the gas flow in a barred potential 
strongly depends on the bar pattern speed. Straight shocks, 
often seen in large-scale bars, occur when the bar is rapidly 
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rotating, i.e. it extends at least to 80% of its corotation. For 
slower bars, the shocks curl around the bar and start form- 
ing a ring (Athanassoula 1992b). Thus in a slowly rotating 
secondary bar, we expect an elliptical flow to develop in- 
stead of straight shocks (Maciejowski 1999). The actual gas 
morphology will obviously be influenced by stars-gas inter- 
action, self-gravity in the gas, and star formation processes, 
none of which is considered in this simple scenario. 



8 SUMMARY 

We have shown that doubly barred galaxy potentials allow 
stable regular orbits that can support the shape of both 
inner and outer bars. These multiply barred galaxies can 
be built as gravitationally self-consistent systems. In dou- 
bly barred galaxies, the most important loop families, which 
serve as their backbones, occupy the same parts of phase 
space as the xi and X2 orbits in the single bar. The loop fam- 
ily corresponding to the xi orbital family in the main bar, 
always remains elongated along the main bar, even though 
the inner loops reside inside the secondary bar. The X2 loop 
family, which at large radii behaves like the X2 orbital family 
of the main bar, follows the secondary bar in its motion at 
smaller radii. No loop family corresponding to the X2 orbital 
family in the secondary bar has been found. 

Despite a limited rmmbcr of models, wc wore able to 
draw conclusions about constraints that the orbital struc- 
ture puts on the shapes, relative sizes and frequencies of the 
bars. The secondary bar can be supported by the X2 loops 
which rotate smoothly with its figure. These loops change 
axial ratio as the bars rotate, and lead or trail the figure of 
the secondary bar, depending on the relative phase of the 
two bars. A self-consistent secondary bar, supported by stars 
trapped around these loops, must pulsate and accelerate as 
it revolves inside the main bar. Since the X2 loops originate 
from the X2 orbits in the potential of a larger bar, the size of 
a self-consistent secondary bar is approximately limited by 
the maximum extent of the X2 orbits along the main bar's 
major axis. Because of those limitations, wo were not able 
to find a self-consistent model of a doubly barred galgixy in 
which the secondary bar remains in resonant coupling with 
the main bar, and extends all the way to its corotation. Our 
nearly self-consistent model contains a slowly rotating sec- 
ondary bar. 

A strong secondary bar can easily disrupt the xi loops 
supporting the main bar; thus there are upper limits for the 
mass and size of the small bar. The inner region of the large 
bar can also be supported by another loop family, which we 
call the bow-tie or br family. The loops supporting one bar 
rotate through the ones supporting the other. Stars can be 
trapped around both these sets of loops, but gas can reside 
on at most one of them. 

The concept of the loop developed in this paper provides 
a powerful tool to study particle orbits in any periodically 
pulsating gravitational potential: in double bars, pulsating 
spheres, and eccentric binaries. It is most promising in the 
last case: the motion of a masslcss test particle under the 
gravitational force of two point masses (the restricted three 
body problem), where only a few isolated closed orbits have 
been found so far. 
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APPENDIX A: LOOPS IN A WEAKLY BARRED POTENTIAL: THE EPICYCLIC SOLUTION 

Here we show how approximations to loops can be found analytically within the epicyclic formalism. We also present the 
solution of a damped epicyclic approximation, which yields excellent approximation to the streamlines of gas flow in a singly- 
barred potential (Lindblad & Lindblad 1994, Wada 1994) . The predictions for gas flow in doubly barred galaxies were explored 
by Maciejewski & Sparke (1997). 

In the epicyclic approximation the motion of a test particle can be decomposed into the 'guiding center' motion, on a circle 
of radius Rq with the angular velocity fl{Ro) corresponding to circular motion in the axisymmetric potential $o, together with 
the epicyclic oscillations resulting from the forcing term $i, and a free oscillation at the local epicyclic frequency k{Ro). In a 
single bar, this formalism gives the approximate orbits that are the counterparts of the xi and X2 families. On these closed 
orbits, the free oscillation is absent. In the same way, the solution of the epicyclic approximation in a double bar can lead 
to finding counterparts of the main loop families. We assume the free oscillation to be absent in the double bar as well. We 
extend the epicyclic solution for orbits in a barred potential (see e.g. Sellwood & Wilkinson 1993, Wada 1994 and Lindblad 
& Lindblad 1994) to double bars. 



Al The epicyclic solution for a Hamiltonian System 

In planar polar coordinates {R, <p) the linearized equation of motion (|l|) yields following equations for _Ri and <pi first order 
corrections to the guiding center motion 

\ dR^ \R„ J dR \Ro,^o 

Rl<pi+2RiRo^lo = , (A2) 

9'P \Ro.<Po 

where flo — Q{Rq). If in addition to the main bar rotating rigidly with angular velocity Qb, there is another, secondary bar 
rotating with angular velocity Qs, the first order term "l>i in the multipole expansion of the potential in the frame of the main 
bar becomes time-dependent, and consists of two bisymmetric terms 

^i{R,(p,t) = *s(_R)cos2v3 + *s(-R)cos(2v?-Ljpt), (A3) 

where cjp — 2(^1$ — f^s), and the time t was chosen such that t = when bars are parallel. The first term on the right-hand 
side corresponds to the main bar, at rest in its own frame, and the second one describes the secondary bar, which rotat es in 



this frame. After substituting Cb,s = ^{d'i! b,s /dr)\Rg and Db,s = b,s / R)\Rq, the right-hand sides of equations (Al) 

and (|A2|), for radial and azimuthal motion respectively, become 

— = Cb cos2((9o + Cs cos(2(y3o — ujpt), (A4) 
OH \Ro,'Po 

1 

--— — — = -Db sm2ipo ~ Ds sm{2ipo - LOpt). (A5) 

In general, we can start a particle orbit at any relative position of the bars, at an arbitrary time ts ^ 0, with its guiding center 
at an arbitrary angle (fios = foit = ts). The angle (fo of the epicycle center at a given time t is then ipo = (fios + (t — is)^ , 
where lj = 2(r2o — ^b)- The two bars induce two characteristic forcing frequencies, so we introduce two time counters, to 
integrate the equations of motion in a way analogous to Wada (1994). The first frequency is lo and the corresponding time t' 
is defined by u jt' — 2ipo. The second frequency is Auj = ui — uip with the correspo ndin g time defined by Acut" = 2ipo — i^pi- 



Equation (A2) for the azimuthal motion, with right-hand side stated explicitly in (A5), can be integrated over time to single 



out ifi and substitute it back to equation (Al) for the radial motion. Thus we have the second order equation for Ri in the 
doubly barred case 

Ri + K Ri = cos Lot + — — cos Atjt , (Ad) 

where tub ~ (jjCb + 2Db^o, rns = A ujCs + 2_Dsf2o, k is the free epicyclic frequency defined by = AD,o{Q,q — A), and A 



is the Oort constant. In equation (A6), one can substitute Ri — Rib + Ris and perform a standard separation of variables. 



which leads to two decoupled equations 

Rig+K^Rig = — ^costjt', (A7) 

Ri^+K^Ri^ = ^cosAtjt", (A8) 
Auj 

of the form exactly like Wada's equation (4), the solution of which we already know. The steady-state solution for Ri in a 
doubly barred potential consists then of two terms and can be written in the form 

= _^l££aB cos[2v3os + io{t — ts)] — ^^^ascos[2v303 + (^{t — ts) — cjpt], (A9) 
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where Aqb.os = ^/Pb,S' <^b,s = PB,s'mB,s, Pb = lu^ — k^, ps ~ Ao;^ — . A similar solution can be written for i^. There 
is an obvious way to extend this method to the arbitrary number of A'^ bars by performing a separation of TV corresponding 

variables. 

Clearly, the orbits in a potential of two independently rotating bars do not close. But from equation ( |A9| ) we can see that 
the time transformation t ^ t + 2tt/uip is equivalent to the initial angle transformation (pos — > fos + iruj/up. This means that 
after time 2tt/u!p, when the two bars return to the same relative orientation, the particle ends up at a place, which at the 
last time was occupied by another particle, with its guiding center at the same radius, but position angle differing by vr^. 
Therefore particles having the same guiding radius interchange positions at consecutive bar alignments, while remaining on the 
curve appropriate to the guiding radius Rq. This curve is the epicyclic approximation to the loop; in this linear approximation, 
a loop is a set of particles which have the same guiding radius, and lack any free epicyclic motion, but respond only to the 
periodic forcing of the two bars. A general orbit with that same guiding center will oscillate about the loop. Loops in the 
epicyclic approximation are presented for Model 1 in Maciejewski & Sparke (1997). 



A2 The case of dissipative motion — the Damped Epicyclic Approximation 

Lindblad & Lindblad (1994) and Wada (1994) showed that if a dissipative term is added to the equations of motion of a 
particle orbiting in a single rotating bar, the pattern of closed orbits in linearized solution closely reflects the gas streamlines 
in flows modeled hydrodynamically. This Damped Epicyclic Approximation gives us a cheap way to preview the hydrodynamic 
simulation. Obviously, this approach cannot reproduce a shock structure, or give the inflow rate. 
Adding a friction term Ff^ic to equation (^ we get 

= -V$ - 2(^3 X r) + \fls\'r + Fmc. (AlO) 

A true friction term Ffric should be proportional to the velocity gradient in the gas, but this information is not available a 
priori, so we must use other forms of the dissipative term. We may postulate that the frictional force is proportional to the 
difference between the particle velocity r and the circular velocity x r at its position in an axisymmetric potential $o, which 
is zero for a particle in circular orbit in that potential. The dissipative force Ffric in the inertial frame can be written as 

Ffric = -2A(r, -fix rO, (All) 

where Vi is the particle's position in the inertial frame, and A is the friction coefficient: the minus sign means that dissipation 
tends to reduce velocity deviations for positive A. The friction described by Lindblad & Lindblad (1994) leads to this same 
expression for Ffric, so we call this form the 'Lindblad friction'. 

The rotation speed Vi in the inertial frame is related to the speed r in a frame rotating with the primary bar by 



Vi — r + J7b X Vi, thus equation (All) becomes 



Ffric = -2A [r - (Jl - JIb) x r] . (A12) 

The bar pattern speed enters the friction formula because the dissipative term depends on both radial and tangential velocity 
deviations from the circular orbit. Wada (1994) used another form of friction, which depends on radial velocity deviations 
only: 

Ffric = -2A(rT)-r/jr|^ (A13) 

In dissipative motion, the two components of particle velocity are coupled, which makes these two forms of friction compatible. 
For Lindblad's friction, we get in a single bar the second order equation for the radial correction Ri 

Ri + 4XRi + 4(^0 + A^ - flQA)Ri = — cos cot' + - smut', (A14) 



where t' is the same as defined for equation (A.6), and n = 2AC When we use Wada's friction defined in equation (A13), the 
second order equation for Ri takes the form 

Ri + 2XRi + 4QoRi{flo - A) = — cosLut'. (A15) 

LU 

The epicyclic frequency k depends on the friction strength for the Lindblad friction, but is the same as in the case with no 
friction when Wada's friction is used. The damping term is half as large for Wada's friction, which damps radial deviations 
only. 



For Lindblad's friction, the dissipative counterpart of the second order differential equation (A6) for Ri in the doubly 
barred case can be written as 

Ri + 4\Ri + K^Ri = — — cos Lot' H — — sin Lut' + — cos Aujt" + sin ALut" , (A16) 

LJ LU AlU AlU 

where the additional coefficients are ris ~ 2XCb and ns = 2\Cs- In presence of friction, we get the following steady-state 
solution for R\ in a doubly barred potential 

-Ri = {ascos[2ipos + Ld{t - ts) — LUpt\+bssm[2ipos + Lu{t — ts) — ijJpt]) 

(as cos[2lpos + Lo{t - ts)] + &s sinpv^os + Lo{t - ts)]) , (A17) 
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where for a dissipative case ub.s = PB,smB,s + <lB,snB,s, bB,s = PB,snB,s — qB,s'mB,s, Qb = 4Acj, and qs = 4AAa;. A 
corresponding solution can be written for ip. Like in the frictionless case, also in the damped epicyclic approximation particles 
with the same guiding radius remain on the same loops. These loops for Model 1 are presented in Maciejewski & Sparke 
(1997). 



APPENDIX B: POTENTIAL OF A PROLATE FERRERS' BAR 



Here we derive the potential of a prolate Ferrers' bar by reducing Pfenniger's (1984) solution for a triaxial bar. It can be 
expressed by a third order polynomial in , , with coefficients Wjm (the reduced version of Pfenniger's coefficients VFjjfc) 
defined by 



u)im{x,y,z) = 



du 



Va2 + m(62 + u) (o^ + uy (62 + u)" 



(Bl) 



The definition of A, and A(A) below, can be found in Pfenniger's paper (see also Binney & Tremaine 1987, pp. 57 and 61). 
Since A = inside the bar, ui's do not depend on location there. Pfenniger's recurrence relations of (All) type reduce to 



COim — (ctJj— 1,T1 

with 



■ i^i,m-i)l{a^ — h^)- 



2i-l^A(A)(a2 + A)(<-i) 
2m^A(A)(62 + A)('"-i) 



- Wi-l,lJ, 



Wl,m,-l], 



Woo 



WlO 



WOl 



In 



/I + sin6'\ 
V cos 61 / 



In 



cos f 
/ 1 + sin e \ 
V cose j 



(a2- 62)3/2' 



tan d sm ( 



■ uJio/2, 



{a? - 62)3/2 

where cos 6 = \J (62 + X)l{o? + A). Then the potential of a prolate Ferrers bar can be written as 

^1 \ poab^ I „ 2 2 2 

2) = ^ — (Woo - 6x 1/ « wi2 

+x^{a:^[3a;2o - x'^'^m\ + 3[y^(2a;ii - j/^wi2 - x^ujix) - 
+2/^{2/^[3wo2 - y^wos] + 3[2;^(2wo2 - 2^a;o3 - J/^woa) - woi]} 
+2^{z^[3a;o2 — a^woa] + 3[a;^(2a;ii — x'u)i\ — z'u)\i) — woi]}) 

The recipe above is equivalent to the formulae given by Perek (1962) and de Vaucouleurs & EYeeman (1972). 



(B2) 

(B3) 
(B4) 
(B5) 
(B6) 
(B7) 



(B8) 
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